#!/bin/csh
# overplot_profiles.chs
###################################################
# Edit 'overplot-profile_files'
# Run the script using: ./1d_overplot.csh field 
# where field can be temperature, velocity, viscosity, composition,
#                    expansivity, conductivity
# NB: the ouput files are stored in the dir1 directory selected in
#     overplot-profiles_files under ../Plots/Comparisons/
###########################################################

#---------------------------
# Options
#---------------------------
set makejpg = 1
set jpgdensity = 200
set deletepdf = 0
set dimx = 10
set dimz = 10

#-------------------------------------------------------------------------
# Read file containing the simulation directory and the profiles to plot
#-------------------------------------------------------------------------
source overplot-profiles_files

#------------------------------
# Set ouput directories
#------------------------------
set dirplots0 = {$dir1}Plots/
set dirplots = {$dir1}Plots/Comparisons/ 

#---------------------------------------------------
# Prepare arrays of profiles and gmt plot options 
#---------------------------------------------------
set prof = ($f1 $f2 $f3 $f4 $f5 $f6 $f7 $f8 $f9 $f10)
set opt = ($opt1 $opt2 $opt3 $opt4 $opt5 $opt6 $opt7 $opt8 $opt9 $opt10)

#---------------------------
# Create output directory
#---------------------------
if (! -d $dirplots0) then
  mkdir $dirplots0
endif

if (! -d $dirplots) then
  mkdir $dirplots
endif

#---------------------------
# Select field to plot
#---------------------------
set field = ${1}
if (${1} !~ temperature && ${1} !~ viscosity && ${1} !~ composition && ${1} !~ conductivity && ${1} !~ expansivity) then
  echo 'The field you want to plot is either not (yet) supported in this routine or has been misspelled'
  exit
endif

#----------------------------
# Set plot origin (in cm)
#----------------------------
set X0 = 3
set Y0 = 5

#--------------------
# Set vertical range
#--------------------
set zmin = 0
set zmax = 1

set output = t.ps

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

####################
# Temperature
####################
if ($field == temperature) then
  echo "Temperature overplot of..."

  set pdfoutput = temperature_profiles_overplot.pdf
  set jpgoutput = temperature_profiles_overplot.jpg

  set title = ""
  set xlabel = "Temperature"
  set zlabel = "Depth"

  set min = 0
  set max = 1
  set delta = 0.2
  set zmin = 0
  set zmax = 1
endif

####################
# Viscosity
####################
if ($field == viscosity) then
  echo "Viscosity overplot of..."

  set pdfoutput = viscosity_profiles_overplot.pdf
  set jpgoutput = viscosity_profiles_overplot.jpg

  set title = ""
  set xlabel = "log@-10@-(Viscosity)"
  set zlabel = "Depth"

  set min = -2
  set max = 4
  set delta = 0.5
  set zmin = 0
  set zmax = 1
endif

####################
# Composition
####################
if ($field == composition) then
  echo "Composition overplot of..." 

  set pdfoutput = composition_profiles_overplot.pdf
  set jpgoutput = composition_profiles_overplot.jpg

  set title = ""
  set xlabel = "Composition"
  set zlabel = "Depth"

  set min = 0
  set max = 1
  set delta = 0.2
  set zmin = 0
  set zmax = 1
endif

####################
# Expansivity
####################
if ($field == expansivity) then
  echo "expansivity overplot of..." 

  set pdfoutput = expansivity_profiles_overplot.pdf
  set jpgoutput = expansivity_profiles_overplot.jpg

  set title = ""
  set xlabel = "Thermal@~\040@~expansivity"
  set zlabel = "Depth"

  set min = 0
  set max = 3
  set delta = 0.5
  set zmin = 0
  set zmax = 1
endif

####################
# Conductivity
####################
if ($field == conductivity) then
  echo "conductivity overplot of..." 

  set pdfoutput = conductivity_profiles_overplot.pdf
  set jpgoutput = conductivity_profiles_overplot.jpg

  set title = ""
  set xlabel = "Thermal@~\040@~conductivity"
  set zlabel = "Depth"

  set min = 0
  set max = 5
  set delta = 1
  set zmin = 0
  set zmax = 1
endif

#------------------------
# GMT plot
#------------------------
psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${delta}:${xlabel}:/a0.2:${zlabel}:WSen:.${title}: -K >! $output
  
set i = 1
while ( $i <= $nprof )
  echo $prof[$i]
  
  set K = -K
  if ($i == $nprof) then
    set K = ""
  endif

  if ($field == temperature) then
    awk 'NR>1 {print $3, 1-$1}' $prof[$i]  >! p
  else if ($field == viscosity) then
    awk 'NR>1 {print log($4)/log(10), 1-$1}' $prof[$i]  >! p
  else if ($field == composition) then
    awk 'NR>1 {print $6, 1-$1}' $prof[$i]  >! p
  else if ($field == expansivity) then
    awk 'NR>1 {print $10, 1-$1}' $prof[$i]  >! p
  else if ($field == conductivity) then
    awk 'NR>1 {print $11, 1-$1}' $prof[$i]  >! p
  endif

  psxy p -R -JX $opt[$i] $K -O >> $output

  @ i = $i + 1
end


#----------------------------------------
# Conversions and cleaning
#----------------------------------------
ps2raster $output -Au -P -Tf -F$pdfoutput
if ($makejpg == 1) then
  ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  mv -f *.jpg $dirplots
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f *.pdf
else 
  mv -f *.pdf $dirplots
endif

rm -f p *.ps 




